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H I I revie^v four decades of numerical simulations of galaxy clusters fo- 

cussing on the attempts to resolve their internal structure. Over- 
merging describes the numerical or physical disruption of dark mat- 
ter halos ^vithin dense environments. This problem was inherent in 
simulations prior to 1998 but can be completely overcome w^ith cur- 
rent algorithms and hardware. We can now resolve many thousands 
of subhalos within clusters and we may have converged on their in- 
ner structure, allow^ing several new tests of the hierarchical structure 
formation model and the nature of dark matter. 



1 Introduction 

The notion that a galaxy cluster forms via gravitational instability within an expanding universe 
was first discussed by van Albada (1961). The first numerical studies of cluster stability with 
10-100 particles were carried out by von Hoerner (1960) and Aarseth (1963). An investigation 
of the virialisation process using several hundred particles was carried out by Henon (1964) and 
Peebles (1970). The assumption being that the "galaxies" had already collapsed prior to the 
turn around of the cluster. From these first simulations it was clear that one could form an 




equilibrium cluster of particles with global properties not too unlike the theorists "ideal" Coma 
cluster. 

White (1976) followed the evolution of 700 unequal mass particles expanding with the Hubble 
flow. In this case it was possible to follow the collapse of groups of particles prior to the formation 
of the cluster. However, these clumps of bound particles were erased during the collapse and 
virialisation of the final system. This work may have been the key motivation behind the 
seminal paper "Core condensation within heavy halos..." (White & Rees 1978). This paper 
argued that cooling and settling of baryons at the centres of dark matter halos was essential 
to the formation and survival of galaxies in dense environments. Although a valid description 
of galaxy formation, it is now apparent that this process is not necessary in order to resolve 
substructure within hierarchical models. 




Figure 1 The hierarchical evolution of a galaxy cluster in a universe dominated by cold dark matter. 
Small fluctuations in the mass distribution are present but barely visible at early epochs. These grow 
by gravitational instability, merging and accretion of mass, eventually collapsing into virialised quasi- 
spherical dark matter halos. This plot shows a time sequence of 6 frames of a region of the universe that 
evolves into a cluster of galaxies. The colours represent the local density of dark matter plotted using a 
logarithmic colour scale. Linear over-densities are darker blue, whereas the non-linear collapsed regions 
attain over-densities of a million times the mean background density and are plotted as yellow/white. 
Each box is 10 Mpc on a side and the final cluster virial radius is 2 Mpc. 

In the 1980's more detailed cosmological models were developed in which the mass of the uni- 
verse was dominated by coUisionless particles. The linear power spectrum of density fluctuations 
was determined by the growth of fluctuations in the expanding universe and the temperature of 
the dark matter particle. The cold dark matter (CDM) model was born (Peebles 1983). 

N-body simulations of cluster formation were used to constrain the hot dark matter model 
(White etal 1984) and most subsequent work has focussed on the cold dark matter model (CDM), 
the currently favoured scenario for structure formation in the universe. Warm dark matter has 
not received a great deal of attention in the literature - primarily because this model could not 
naturally produce a closure density of dark matter. Attention has now shifted from standard 
CDM models, which fail for a variety of reasons, to a Lambda dominated CDM universe. Since 
the universe is clearly not as simple as one might wish, alternative models such as warm or self 
interacting dark matter may become more popular. Perhaps the widespread acceptance of a 
positive Lambda term stems mainly from the fact that the most hierarchical model fails without 
its inclusion. Although a hierarchical universe can explain many fundamental observations, from 
the microwave background to clustering at high redshifts, we still know next to nothing about 
the physical nature of the dark matter. 

Many of the numerical simulations during the 1980's focussed on a comparison between 
models and observations of large scale structure. Studies of individual halo properties were 
taken from large cosmological volume simulations. For example, Frenk etal (1985) and Quinn 
etal (1986) analysed the structure of dark matter halos that contained just a few thousand 
particles per object. Code development, "volume renormalisation" and faster computers lead 
to the ability to simulate individual clusters with more than 10^ particles and with a force 
resolution of approximately 1-2% of the virial radius (Dubinski & Carlberg 1991, Warren etal 
1992, Summers etal 1996, Tormen etal 1996, 1997, 1998). These simulations still produced 
smooth dark matter halos with very little surviving dark matter subhalos - although interesting 
results were found concerning the inner structure of halos and the formation and accretion 
history of halos. 

Warren etal and Carlberg were the first to resolve the inner 10% of a halos virial radius 
and claimed a power law slope p(r) oc r~^. Later simulations by Navarro etal (1997) revealed 
the remarkable scaling properties of CDM halos across a wide range of mass scales, from dwarf 
galaxies to galaxy clusters. Determining the central cuspy profile of dark matter halos is an 
important but difficult computational problem and we have seen recent progress in this area 
from many groups. 

Analytic work concerning the overmerging problem began with White & Rees (1978) who 
discussed the various numerical and physical processes that could lead to the loss of substructure. 
Carlberg (1994) argued that particle-halo heating was the cause of overmerging, a result also 
claimed by van Kampen (1995). However the timescale for this process was at least a Hubble 
time whereas overmerging appeared on an orbital timescale. Moore etal (1996) argued that the 
resolution at that time was sufficient to overcome relaxation and tidally accelerated disruption 
by the finite mass background particles. 

The physical processes of tidal disruption via halo-halo collisions and by the global cluster 



potential take place on a more rapid timescale and could explain the loss of substructure in the 
simulations prior to 1998 (Moore etal 1996 - the investigation of the importance of this process 
lead to the notion of galaxy transformation in clusters by tidal heating - galaxy harassment.) 
Tidal disruption is rapidly enhanced if halos are poorly resolved in their central regions (c/ 
Figure 2). The survival of dark matter substructure depends critically on the central density 
structure which in turn depends on the force and mass resolution in a simulation. For example, 
if halos have resolved singular density profiles p{r) oc r~^, then it is physically impossible to 
entirely disrupt a satellite halo. Only through physical merging can such halos be lost. 

The key motivation behind the Seattle HPCC group ( http://www-hpcc. astro. Washington. edu/| ) 



was to construct a high performance parallel treecode (PKDGRAV) that could simulate the for- 
mation of cosmic structures with high resolution. The aim was to resolve suh — galactic halos 
with thousands of particles and ~ kpc force resolution. After a long testing program a set of 
accuracy parameters were chosen that allowed the correct growth of fluctuations to be followed 
on large scales, as well as accurately following the orbits of particles in the most dense regions. 




Figure 2 The same cluster siniulated in Figure 1 but with two different values of the softening length 
and keeping the particle mass fixed. The left panel shows a close up view of the inner 500 kpc of the 
last frame of Figure 1. In this case the softening was 0.2% of the virial radius. The right panel shows 
the same region of the same cluster but simulated with a softening length of 1.5% of the virial radius. 
The lack of substructure halos in the right panel demonstrates that softened halos are easily disrupted 
by tidal forces. 

The first simulation to finish was the 1.3 million particle "Virgo" cluster simulation presented 
in Moore etal (1998) and analysed by Ghigna etal (1998). This was a simulation of the formation 
of a single galaxy cluster within a standard CDM universe (Figure 1). The box size was 100 
Mpc but a series of nested shells of different mass particles allowed the cluster to be resolved 
at high resolution. Many hundreds of dark matter subhalos could be found orbiting within 
the cluster - the overmerging problem was resolved. A higher resolution study of the same 
cluster demonstrates that the central density profile appears to have converged as have the 
global properties of the subhalos within the cluster (Ghigna etal 2000). (Similar high resolution 
studies were being carried out at the same time by several other groups (c.f. Klypin etal 1999, 
Cohn 1999, Okamoto & Habe 1999, Jing 2000, Fukushige & Makino 2000). 

It is interesting to see how numerical calculations are scaling with respect to algorithmic and 
computational developments. The cluster simulation by Peebles (1970) contained 300 particles 



and had a force resolution of order 100 kpc. The highest resolution simulation published to 
date is by Ghigna etal (2000) in which a cluster halo was simulated with 10 million dark matter 
particles and force softening of 1 kpc. Roughly a factor of 10^ times the computational cost in 
a timescale of 30 years - an increase in speed that arises from both algorithmic and hardware 
developments. Four major development have lead to this performance increase: grid based codes 
or treecodes reduce the work of long range force calculations, multistepping saves moving all the 
particles on the same short timestep as the handful of particles in dense regions, faster cpu's, 
parallel codes can scale up to 90% efficiency on large numbers of nodes. Figure 3 shows the "N" 
versus year plot for several cluster simulations over the last 38 years. 
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Figure 3 A selection of dark matter only cluster simulations. I have attempted to plot the number of 
particles within the final virial radius against the publication year. The dashed line illustrates "Moore's 
law" (named after Gordon Moore who founded Intel and noted the rapid increase in the density of 
transistors on a silicon chip with time) which roughly equates to a doubling of cpu speed every three 
years. This type of plot suffers from many biases and is not meant to be taken seriously. The fact that 
recent simulations lie on this curve does not imply that algorithmic advances are not important quite 
the opposite since the force softening is reduced typically in proportion to N~^'^ which requires more 
computational work. 



We now have a host of new observational tests of the hierarchical structure formation 
paradigm that include: the extent of galactic halos in clusters or satellites in galactic halos, 
the orbital distribution of cluster /satellite galaxies, the number of satellites as a function of 
their circular velocity or mass, the spatial and velocity distribution of satellites, the central 
density profiles of clusters and galaxies are close to r~^'^. Unfortunately, observational tests of 
these model predictions on cluster scales are non trivial and just a few exist in the literature 
(Tyson etal 1998, Natarajan etal 1999, Smith etal 2000). 
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